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Abstract 


All  commonly  used  failure  models,  both  univariate  and  multivariate,  are  indexed 
by  a  single  scale,  typically  time.  This  thesis  considers  the  case  of  indexing  a 
failure  model  by  two  scales,  namely  time  and  usage.  A  family  of  models  that  is 
indexed  by  these  two  scales  of  measurement,  and  attempts  to  capture  the 
dependence  between  them,  is  developed  and  an  example  given  in  detail. 
Extensions  to  this  basic  model,  and  problems  associated  with  them,  are  explored. 
As  an  aside,  an  idealised  stochastic  process  for  modelling  usage  is  introduced. 
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Failure  models  indexed  by  time  and  usage 
by  Simon  P.  Wilson 


1  Introduction 

A  glance  at  any  book  or  journal  in  the  field  of  reliability  and  survival  analysis 
will  show  that  much  time  and  effort  has  been  devoted  to  probabilistic  models  for 
failure  of  an  item  or  system  of  items.  One  can  choose  to  model  the  time  to 
failure  of  an  object  from  a  wide  variety  of  univariate  probability  distributions 
such  as  the  exponential,  Weibull  or  gamma.  The  times  to  failure  of  a  system  of 
units  can  be  modelled  by  one  of  many  multivariate  distributions;  see  Hougaard 
(1987)  for  a  review  of  some  of  the  more  commonly  used.  All  these  distributions 
have  a  feature  in  common,  namely  that  they  are  indexed  on  a  single  variable  or 
scale  of  measurement.  Often  this  variable  is  time  but  may  be  number  of  cycles, 
amount  of  exposure  to  a  hazardous  material;  any  variable  that  the  modeller 

thinks  relevant  in  fact. 

Whilst  these  models  are  widely  used  and  axe  very  versatile,  there  are  certain 
situations  where  one  would  like  failure  of  an  item  to  be  indexed  by  more  than 
one  variable.  In  this  case  the  existing  models  are  not  satisfactory.  A  classic 
example  is  in  the  area  of  automobile  warranties,  where  a  manufacturer  will  agree 
to  repair  or  replace  a  vehicle  if  it  fails  within,  say,  3  years  or  30,000  miles  from 
production.  If  a  manufacturer  wishes  to  investigate  the  expected  number  of 
warranty  claims  or  tries  to  find  a  warranty  that  maximises  his  profit  margin, 
then  it  is  necessary  to  have  a  failure  model  for  automobiles  that  is  in  terms  of 
time  to  failure  and  number  of  miles  to  failure.  Another  example  would  be 
measuring  the  number  of  deaths  amongst  mineworkers  in  terms  of  age  at  death 
and  amount  of  exposure  to  hazardous  dust  whilst  at  work.  In  both  these 
examples  the  two  scales  of  measurement  are  obviously  dependent  -  what  is 
needed  is  a  failure  model  indexed  jointly  by  the  two  scales.  In  this  thesis  we 
assume  these  two  scales  are  time  and  usage,  but  this  work  could  equally  be 
applied  to  time  and  some  other  index  like  exposure  to  hazardous  dust. 
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To  the  best  of  my  knowledge  almost  no  work  has  been  published  that  looks  at 
this  situation.  That  failure  can  be  represented  by  more  than  one  scale  of 
measurement  has  been  recognised,  see  Farewell  and  Cox  (19/9)  or  Oakes  (1988). 
In  these  papers  the  authors  consider  the  problem  of  reducing  several  scales  of 
measurement  down  to  one  that  is  most  informative  in  some  sense,  with  the 
starting  scales  becoming  time-dependent  covariates. 

The  innovative  nature  of  the  idea  here  must  be  emphasised;  this  topic  has 
received  almost  no  attention  in  the  literature  despite  the  fact  that  such  models 
arise  naturally  in  many  situations  of  interest.  A  naive  solution  to  the  problem  of 
having  a  suitable  model  would  be  to  select  one  of  the  available  bivariate  failure 
models  and  then  attempt  to  justify  its  use  as  a  time/use  failure  model  in  some 
manner.  However,  this  is  not  the  approach  taken  in  this  thesis;  instead, 
assumptions  are  made  concerning  the  dependence  between  the  two  scales  of 
measurement  and  a  model  is  derived  from  those  assumptions.  It  is  hoped  that 
the  result  is  a  model  with  more  justification  than  one  developed  from  the  former 
approach,  which  would  employ  a  model  that  was  really  designed  to  describe  the 
joint  failure  of  two  units  in  terms  of  one  scale  rather  than  the  failure  of  one  unit 
in  terms  of  two  scales. 

This  thesis  is  split  into  5  further  sections.  Section  2  introduces  a  family  of  models 
indexed  by  time  and  usage  then  section  3  considers  an  example  of  this  model  in 
some  detail.  Section  4  considers  generalisations  of  the  model  and  some  problems 
that  may  be  encountered  with  them,  whilst  section  5  looks  at  a  more  realistic 
stochastic  process  to  model  usage.  Section  6  concludes  the  thesis  and  outlines 
possible  areas  for  future  research  in  this  topic. 
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2  A  bivariate  failure  model  indexed  by  time  and  usage. 

The  aim  of  this  section  is  to  introduce  a  generic  bivariate  model  for  failure  in 
terms  of  both  time  and  usage. 


Consider  some  unit  or  object  like  a  car,  a  light  bulb  or  any  item  which  can 
accumulate  use;  in  the  case  of  an  automobile  this  use  will  be  mileage,  in  the  case 
of  a  light  bulb  it  may  be  number  of  times  switched  on  or  the  amount  of  time 
that  it  is  lit.  To  begin  with,  we  observe  that  this  unit  may  age  over  time,  even 
though  it  is  not  used,  due  to  rusting,  humidity,  adverse  temperature  and  other 
environmental  effects.  These  effects  are  dynamic,  that  is  change  over  time. 
Singpurwalla  and  Youngren  (1990)  have  shown  that  if  w(t)  denotes  the  failure 
rate  of  the  unit  under  ideal  (or  laboratory)  storage  conditions,  and  if  the 
dynamics  of  the  storage  environment  are  described  by  a  stochastic  process  known 
as  a  “gamma  process”  with  parameters  a(t)  and  b,  then  the  failure  rate  of  the 

unit  under  storage  is 

(2.1)  r(t)  =  a(0  ln  (  1  +  ) 

If  _  a  an(i  w(t)  =  u  then  r(t)  is  a  constant  and  so  ageing  under  storage 
would  be  that  provided  by  an  exponential  distribution.  Similar  results  are 
obtained  when  the  gamma  process  is  replaced  by  a  “shot-noise  process  . 

From  (2.1)  we  have  a  failure  rate  for  the  unit  when  it  is  not  used.  Of  course  one 
can  choose  any  failure  rate  under  non-use  that  one  likes;  the  point  of  (2.1)  is  to 
demonstrate  that  models  have  been  developed  in  the  literature  to  meet  this 
requirement.  Now  let  us  put  the  unit  to  use,  and  denote  the  amount  of  usage  at 
time  t  as  M(t).  We  use  the  notion  of  proportional  hazard  (Cox  (1972))  and 
suppose  that  the  effect  of  usage  is  to  increase  the  conditional  failure  rate 
r(t|M(t))  of  the  unit  from  r(t)  to  say  r(t)  +  /?M(t),  where  /?  >  0  is  a  constant; 
other  forms  are  possible.  We  are  uncertain  about  the  usage  at  time  t  and  so  M(t) 
is  assumed  to  be  some  positive  and  non-decreasing  stochastic  process,  in 
particular  a  Poisson  process  with  intensity  function  A(t).  Intuitively  A(t) 
measures  the  rate  of  use  of  the  item  over  time,  and  figure  1  gives  three  possible 
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functions  that  are  candidates  for  A(t):  respectively  a  constant  rate  of  use,  cyclical 
use  and  a  rate  of  use  that  decreases  as  the  unit  becomes  old.  Observe  that  this 
makes  M(t)  a  discrete  quantity  taking  values  in  the  non-negative  integers, 
although  one  is  at  liberty  to  consider  other,  possibly  real-valued,  non-decreasing 

stochastic  processes  to  model  M(t). 


Figure  1  Three  candidate  functions  for  A(t) 

With  these  assumptions  we  have  tried  to  model  the  physics  of  the  relationship 
between  time  and  usage  to  failure  in  a  simple  but  intuitive  way.  In  any  event,  we 
are  now  able  to  derive  a  bivariate  density  for  failure  of  the  unit  in  terms  of  time 
and  use  based  on  these  assumptions.  The  derivation  is  quite  straightforward,  and 
will  only  be  outlined  below;  readers  interested  in  the  full  technical  details  of  the 
derivation  are  referred  to  the  appendix. 
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Let  T  denote  time  to  failure  and  U  usage  to  failure.  The  first  step  is  to  reason  as 
follows 


(2.2)  f'p  jj(t,  u)  =  density  of  failure  at  time  t  and  use  u 

=  density  of  failure  at  time  t  and  with  use  at  time  t  of  u 

=  fT,M(t)^’  u) 

=  fT|M(t)^  u)  fM(t)(u) 

_  f  ft,  u)  e-A^) 

“  fT|M(t)(t>u)  u! 

where  A(t)  =  }  A(s)  ds,  since  M(t)  is  a  Poisson  process. 

0 

Using  the  well  known  identity  connecting  any  continuous  probability  density  f(t) 
with  its  failure  rate  h(t) 

t 

f(t)  =  h(t)  exp  (  -  j  h(s)  ds  ) 


we  can  say 


t 

(2.3)  fX|M(t)(t|  M(t)=u)  =  (r(t)+/?u)  x  exp  (-  j  (r(s)+/?M(s))  ds  |  M(t)=u  ) 


where  Jf(s)  ds  |  f(t)=u  denotes  the  integral  of  f(s)  given  that  f(t)-u. 


By  using  this  expression,  we  can  follow  the  method  outlined  by  Lemoine  and 
Wenocur  (1986)  and  arrive  at 


,  .  .  ,  (r(t)  +  j8u)  e 

(2.4)  u*  “  (  A(t)  )u 


■R(t) 


A(s)  e 


-/?(t-s) 


ds  )U 


where  R(t)  =  J  r(s)  ds,  the  cumulative  hazard. 

0 
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Now  incorporate  (2.4)  into  (2.2)  to  obtain 


(2.5)  u) 


(r(t)  +  pn)  e 


,-R(t)  -A(t)  ,  t 


(  |  A(s)  ds  )  u 


U! 


'  where  t  >  0,  u  =  0,  1,  2,  ... 


It  must  be  emphasised  that  this  is  not  a  bivariate  density  for  failure  in  the  usual 
sense  of  modelling  the  joint  failure  of  two  items,  but  rather  describes  the  failure 
of  a  single  item  in  terms  of  time  and  use.  There  is  flexibility  m  this  model,  as 
one  is  at  liberty  to  assign  any  intensity  of  use  function  A(t)  and  failure  rate  for 
storage  r(t)  that  seems  plausible,  and  indeed  to  alter  the  relationship  r(t)+/?M(t) 
between  them.  Note  that  the  integral  appearing  in  (2.5)  is  the  convolution  of  A(s) 
with  e^s,  so  can  be  found  via  Laplace  transforms  if  direct  integration  proves 

troublesome. 


The  above  density  enjoys  some  analytic  tractability,  even  in  the  most  general 
case.  For  instance,  the  marginal  of  T  is 


(2.6)  fT(t)  =  (r(t)+/? 


A(s)  ^  ds)  x  exp  {-R(t)  - 


A(s)  (l  -  e^^)  ds} 


t  >  0. 


This  gives  a  marginal  survival  of 

(2.7)  F(t)  =  P(T>t)  =  exp  {  -R(t)  -  j  A(s)  (l  -  )  ds  }  t  >  0. 

0 

It  is  then  an  easy  matter  to  find  the  marginal  failure  rate  as 

(2.8)  rT(t)  =  ^  =  e_',<t'S)  ds  1 2  °  ' 

Equations  (2.6),  (2.7)  and  (2.8)  defines  a  probability  for  time  to  failure  that 
considers  both  environmental  and  usage  stresses  on  a  unit,  and  as  such  is  of 
merit  by  itself  as  a  univariate  model. 
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Given  (2.5)  and  (2.6)  we  have  the  conditional  density  of  usage  to  failure  given 
time  to  failure  as 


(2.9)  P(U=u|T=t) 


r(t)  +  flu  {v{t))n  -r}( t)  u  - 

r(t)  +  0t)(  t)  u! 

where  rj(t)  =  J  A(s)  e  ^  ^  ds 

0 


Note  that  this  is  a  different  density  from  M(t),  usage  at  time  t,  unless  (3-0. 
Whereas  the  expected  value  of  M(t)  is  A(t),  the  expected  value  of  U  given  T=t 
can  be  shown  to  be 


(2.10)  E(U  I  T=t  )  =  ,(t)  x  (  1  +  r(t)  l  ^(t )  )  • 

which  is  seen  to  be  equal  to  A(t)  if  and  only  if  (3=0.  In  general,  the  marginal  of  U 
is  difficult  to  calculate,  although  as  we  see  in  Section  3  it  can  often  be  obtained 
in  special  cases. 


Uncertainty  about  the  parameter  (3  might  be  described  by  placing  a  prior 
distribution  on  it.  Uncertainty  about  the  functions  r(t)  and  A(t)  could  be 
described  by  making  them  stochastic  processes;  thus  M(t)  would  become  a 
doubly  stochastic  or  conditional  Poisson  process.  Section  3  looks  at  these  ideas 
for  a  particular  example  of  the  model. 
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3  Example 


In  this  example  we  consider  the  simple  situation  of  the  failure  of  an  electrical 
switch  in  terms  of  T,  the  time  to  failure  and  U,  the  number  of  times  the  switch 
is  used.  We  assume  that  under  storage  conditions  the  switch  ages  according  to  an 
exponential  distribution  with  rate  a  (r(t)  =  a  in  the  notation  of  section  2).  Each 
time  the  switch  is  used  this  failure  rate  is  incremented  by  an  amount  0  because 
of  the  extra  stress  induced  upon  it  (so  r(t|  M(t))  =  a+/?M(t)  in  the  notation  of 
section  2).  Finally,  we  assume  that  use  of  the  switch  occurs  as  a  homogeneous 
Poisson  process  (that  is  to  say  the  time  between  uses  of  the  switch  are 
independent  and  identically  distributed  exponential  random  variables)  with  rate 
A.  Under  these  assumptions  we  may  apply  our  model  (2.5)  and  obtain  the  density 
for  failure  at  time  t  and  number  of  uses  u  as 


(3-D  fT,p(.,u)  =  (“+<’U);,'(°—  (^(l-e^))U;  t>0,  u  =  0,  1,  2,  ... 


The  density  is  continuous  in  t  but  only  defined  for  u  an  integer,  so  consists  of  a 
sequence  of  univariate  functions,  continuous  in  t,  equally  spaced  along  the  U-axis 
at  u=0,  1,  2,  ...  .  The  plots  of  such  a  density  will  be  referred  to  as  bivariate 

density  sheets. 


The  formulae  given  by  (2.6)  to  (2.10)  can  be  applied  by  letting  r(t)  =  a  and  A(t) 
—  for  example  the  marginal  survival  to  time  t,  and  marginal  failure  rate,  of 

the  switch  axe 


(3.2) 


Frp(t)  =  exp  {  -(a  +  A)t  +  ^(1  -  e  /?t)} 

Trp(t)  =  a  +  A  (1  -  e  ,  t  >  0  . 


In  this  special  case  we  can  exploit  the  homogeneity  of  the  Poisson  process  to 
obtain  the  marginal  for  use  to  failure: 


(3.3) 


P(U  =  u)  = 


(a  -f  ftu)  Au 
H  (A  +  q  +  /?i) 

i=0 


u 


0,  1,2,  ... 
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We  also  have  an  expression  for  the  distribution  function. 

„  An(a+0n)  f  ^  i  ,  p  (  n 

(3.4)  P(T  <  t,  U  <  u)  =  pa.  ni  |  .E  ^  l  i  ) 


i  (  n  ^  0_Z  exp{-(A+a+0i)t}) 
A+a+0i 


Figure  2  is  a  plot  of  the  density  (3.1)  when  a=0.2,  0-0.2  and  A-1.5. 

The  marginal  moments  for  T  and  U  can  also  be  written  down  in  this  case 
although  only  rather  messily  in  terms  of  the  incomplete  gamma  function  7(a,  n). 


Let  x  =  A//J  and  y  =  o/d-  Then 


(3.4) 


E(T)  =  ^  x-(x+y)  X  y(x+y,  X) 


and 


(3.5)  E(U)  =  exp(x)  x^X+y^  x  (  7(2+x+y,  x)  -  (l+2x+y)  7(l+x+y,  x)  + 

x(l+x+y)  7(x+y,  x)  )  . 

Now  we  consider  the  placing  of  priors  on  the  parameters  (a,  0,  A)  of  (3.1),  and 
address  the  question  of  using  data  to  update  our  prior  distributions.  The  first 
observation  is  that  if  a  and  0  are  independent  of  A  then  the  gamma  distribution 
is  a  conjugate  prior  for  A.  This  is  easy  to  see  from  (3.1)  by  noting  that 


L(A  |  t,  u)  =  f(t,  u  |  A) 
oc  Au  e'tA 

Suppose  A  has  a  gamma  distribution  with  shape  ^  and  scale  *r  Given  n 
(independent)  observed  time-use  failure  pairs  (t1?  u^,  (tn,  un)  we  have  a 
posterior  for  A  : 

(3.5)  *(A  |  (t1?  ux),  (tn,  un))  oc  tt(A)xL(A  |  (tx,  uL),  (tn,  un)) 
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so  that  the  posterior  of  A  is  gamma  with  shape  i^+Euj  and  scale  ^E^- 

The  likelihoods  for  a  and  0  do  not  appear  to  possess  a  recognisable  conjugate 
prior.  However,  in  the  case  of  a  a  posterior  can  still  be  written  down  in  closed 
form,  albeit  an  increasingly  complex  one  as  we  increase  the  amount  of  data,  if  we 
assume  a  gamma  prior  on  a.  Assume  a  is  independent  A.  The  likelihood  of  a  can 

be  written  : 


L(a  |  0,  t,  u)  a  (a+/?u)  e'at 

Assume  that  the  prior  on  a  is  a  gamma  with  shape  and  scale  $2.  Then  the 
posterior  of  a  is  of  the  form 

(3.6)  x(o  |  0,  (t1,u1),...,(tn,un))  «  »(a)xL(tt  |  0,  (t1,u1),...,(tn,un)) 

«  e(#2+S‘i)“x  a”2'1  ft  («+/*;) 

i=l 

The  normalising  constant  here  will  be 

(3.7)  j  e'(#2+  E  ‘i*  “  x  /2'1  ft  («+<M  d“ 

o  1=1 

Expanding  the  product,  we  see  that  the  integrand  is  a  finite  sum  of  powers  of  a 
multiplied  by  exp{-(<?2+ 50  (3J)’  the  normalising  constant  for  the 
posterior  of  a  given  by  (3.6),  is  just  a  finite  sum  of  gamma  functions.  Thus  we 
can  give  the  exact  posterior  of  a.  This  posterior  is  dependent  upon  0 ,  so  even  if 
we  assume  a  and  0  to  be  a  priori  independent,  they  are  dependent  in  the 

posterior. 

Similar  calculations  for  0  are  less  straightforward,  due  to  the  more  complex 
nature  of  the  likelihood.  There  is  no  closed  form  posterior  for  0  for  any  of  the 
commonly  used  priors,  and  it  appears  that  one  would  have  to  employ  simulation 
or  some  other  approximation  to  obtain  a  posterior.  If  0  were  known  or  assumed 


10 


as  some  value,  then  (3.5)  and  (3.6)  show  that  one  could  proceed  to  calculate 
exact  posteriors  for  »  and  A  given  that  they  had  independent  gamma  priors.  The 
independence  of  o  and  A  would  remain,  so  these  calculations  could  be  done 
separately.  If  some  prior  was  placed  on  0,  then  we  would  need  to  calculate  a 
joint  posterior  for  .  and  (3  because  of  the  dependence  that  is  introduced  between 

them. 
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4  Using  other  stochastic  processes  to  model  usage 


There  are  several  ways  by  which  the  model  outlined  in  Section  2  can  be 
generalised.  One  that  comes  straight  to  mind  is  to  model  cumulative  usage  as 
some  stochastic  process  other  than  Poisson,  perhaps  some  process  taking  values 
in  [0,oo)  instead  of  the  non-negative  integers.  There  are  numerous  stochastic 
processes  that  one  could  look  at,  and  intuitively  any  process  M(t)  that  satisfies 

i)  M(0)  =  0 

ii)  M(t)  >0  Vt 

iii)  M(t)  is  non- decreasing 

is  a  possible  candidate  to  model  usage. 

Given  that  we  decide  on  some  process  M(t),  then  from  (2.2)  and  (2.3)  we  see 
that  to  obtain  a  density  for  failure  in  time  and  usage  requires  us  to  know  two 
things  about  M(t) 


a)  The  density  of  M(t)  and 


Recall  that  we  assume  the  failure  rate  of  a  unit  when  it  is  not  used  is  r(t).  The 
effect  of  usage  m  is  to  increase  the  failure  rate  to  r(t)+/?m,  where  0  is  a  constant. 
In  other  words  : 

r(t  |  m)  =  r(t)+/?m 

Since  usage  is  assumed  to  be  a  stochastic  process  M(t),  we  rewrite  the 
conditional  failure  rate  as: 

r(t|  M(t))  =  r(t)+/?M(t) 

What  we  desire  is  the  density  of  failure  at  time  t  and  usage  u.  Let  T  be  the 
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random  variable  denoting  time  to  failure,  and  let  U  be  the  random  variable 
denoting  usage  to  failure.  As  in  section  2,  we  reason  as  follows: 

(4.1)  f(t,  u)  =  density  of  failure  at  time  t  and  usage  u 

=  density  of  failure  at  time  t  and  with  usage  at  time  t  of  u 

=  f(t,  M(t)=u) 

=  f(M(t))  x  f(t  |M(t)) 


Using  our  conditional  failure  rate  r(t|  M(t))  we  can  say 

t 

(4.2)  f(tj  M(t)=u)  =  (r(t)-h/?u)  x  exp  (  -  J  r(t)+£M(s)  ds  |  M(t)=u  ) 

0  t 

=  (r(t)+/?u)  e'R^  x  exp  (-  /?  j  M(s)  ds  |  M(t)=u  ) 


where  R(t)  =  j  r(s)  ds 

0 

Combining  equations  (4.1)  and  (4.2)  we  see  that  the  required  joint  density  is  of 
the  form 

t 

(4.3)  f(t,  u)  =  f(M(t)=u)  x  (r(t)+/?u)  eR^  exp  (-/?  j  M(s)  ds  |  M(t)=u  )  . 


(4.3)  provides  a  most  general  expression  for  a  class  of  densities  indexed  by  time 
and  usage.  By  making  M(t)  a  Poisson  process,  we  obtain  the  model  of  section  2. 


What  about  M(t)  as  another  stochastic  process?  I  have  studied,  as  a  first 
example,  the  case  of  modelling  M(t)  as  a  homogeneous  compound  Poisson 
process.  This  can  be  described  as  a  homogeneous  Poisson  process  where  the  jump 
at  each  event  is  no  longer  1,  but  some  random  variable.  I  considered  the  case 
where  the  jump  sizes  are  independent  and  identically  distributed  exponential 
random  variables.  In  other  words,  if  N(t)  denotes  the  number  of  events  at  time  t 
(which  is  a  homogeneous  Poisson  process)  and  Ui  are  lid  exponential  then 

0,  if  N(t)=0 

N(t) 

£  U  if  N(t)  >  1 

i=l 
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This  process  satisfies  the  three  criteria  i),  ii)  and  hi).  A  typical  realisation  of  this 
process  is  given  by  figure  3.  M(t)  can  now  take  any  non-negative  value.  To  give 
this  model  some  practical  explanation,  imagine  that  M(t)  is  the  number  of  miles 


t 

Figure  3  Typical  realisation  of  a  compound  Poisson  process  M(t) 

an  automobile  has  driven.  Under  this  model  we  are  then  assuming  that  the  car  is 
used  according  to  a  homogeneous  Poisson  process,  and  that  each  time  it  is  used 
the  number  of  miles  driven  is  an  exponential  random  variable.  Amongst  other 
things,  this  means  that  we  assume  each  time  between  uses  of  the  vehicle  is 
independent.  Since  this  process  is  a  ’jump’  process,  we  are  also  tacitly  assuming 
that  these  miles  are  accumulated  instantly,  whereas  of  course  in  reality  they  are 
accumulated  over  some  period  of  time.  If  we  assume,  as  we  did  in  the  example  of 
Section  3,  that  the  conditional  failure  rate  is  a  +  /?M(t),  then  we  are  ready  to 
derive  the  joint  density  via  the  same  method  employed  in  Section  2.  M(t)  is 
compound  Poisson;  let  the  usage  events  occur  as  a  Poisson  process  with  rate  A 
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and  let  each  jump  in  the  process  be  exponential  with  mean  1.  The  density  of 
M(t)  can  be  computed;  it  is  continuous  on  (0,oo)  and  has  positive  probability 
mass  at  0,  since  there  is  always  a  positive  probability  that  N(t)=0. 


The  density  of  M(t)  is  given  by: 


(4.5)  P(M(t)  =  0)  =  e 


-At 


5  fM(t)(u)  =  e 


-U-At  f*  M*  u°'1 

n!  (n- 

n=l 


(n-1)! 


u  >  0 


^(t|  u).  Observe  from  (2.3)  that  this  is 

t 


It  remains  to  calculate  ^'j’|M(t)'' 

(t|  u)  =  (a+/?u)  e'at  x  exp  {  -  /?  [  M(s)  ds  |  M(t)=u  } 


(4-6)  fT|M(t)^ 


0 


To  find  Jq4  M(s)  ds,  the  area  under  M(s)  between  0  and  t,  given  that  M(t)  u 
was  quite  straightforward  in  section  2  when  M(s)  was  Poisson  but  is  now  much 
more  difficult.  I  have  been  unable  to  obtain  a  closed  form  for  (4.3)  when  M(t)  is 
compound  Poisson,  but  by  conditioning  on  N(t)  and  combining  with  (4.5)  I  can 
say  that  for  t  >  0  and  u>0  there  is  a  joint  continuous  density  of  the  form 


(4.7) 


fTjU(t,u)  =  («  +  0u) 


X 


00 


Y  An  In(t,  u;  0) 

n=l 


where  In(t,  u;  0)  is  a  2n-dimensional  integral  given  by 


in(45 u;  P) 


exp  {  0  Y  si  ui  }  d(?’ 

i=l 


and  =  {  (s,u)  e  R2n  |  Uj  >  0 


£  u-  =  u  ;  0  <  s1  <  ...  <  sn  <  t  } 
i=l 


As  with  (2.5),  the  details  of  the  derivation  of  (4.7)  are  rather  cumbersome  but 
present  few  difficulties  and  have  been  relegated  to  the  appendix.  This  integral 
has  finite  value,  and  of  course  can  be  evaluated  numerically.  It  is  possible  to 
analytically  integrate  out  (Sl,...,sn),  leaving  an  n-dimensional  integral  over 
(ulr..,un)  to  calculate  numerically.  We  are  integrating  over  the  simplex  £^=11, 
and  so  can  reduce  this  integral  to  an  (n-l)-dimensiona 1  problem.  This  is  still  a 
rather  forbidding  task  for  even  moderate  values  of  n,  and  to  become  practicable 
further  simplification  of  (4.4)  will  be  necessary. 
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As  there  is  positive  probability  mass  at  M(t)=0,  so  there  is  mass  at  U=0  (that  is, 
along  the  T-axis)  in  the  joint  density.  We  can  show  that 


(4.8)  P(U=0)  =  ■ 

It  is  distributed  along  this  axis  with  mass  function 

fTU(t,  0)  =  «  e-(«+A>*  t>0 

Although  the  joint  density  appears  to  be  beyond  a  closed  form  representation,  by 
following  the  method  of  Lemoine  and  Wenocur  (1986)  one  can  obtain  the 
marginal  survival  for  time  to  failure  as  : 

A 

(4.9)  P(T>t)  =  (1+Z^t)^  e'^+A)t  t>0 


The  marginal  failure  rate  is 

(4.10)  rT(t)  =  a  +  A  t~° 

We  notice  a  similarity  between  (4.10)  and  the  marginal  failure  rate  of  the  model 
example  in  section  3,  given  by  (3.2).  Recall  that  the  model  there  simply  assumed 
M(t)  as  a  homogeneous  Poisson  process  with  intensity  A.  In  both  cases  we  have  : 

rT(0)  =  a  ;  hm  rT(t)  =  a+A  ;  rT(t)  is  strictly  increasing. 

The  density  specified  in  (4.7)  is  more  easily  approximated  via  simulation.  The 
integral  in  (4.6)  can  be  easily  simulated  by  computer  and  this  enables  us  to 
approximate  a  density  value  for  any  given  t  and  u.  If  no  further  simplification  of 
(4.7)  is  possible,  then  numerical  simulation  is  the  easier  and  quicker  option  over 
numerical  integration.  Details  of  a  simulation  plan  are  to  be  found  in  the 
appendix.  Figure  4  overpage  is  a  plot  of  the  simulated  joint  density  for 

a=/?=A=l. 
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In  conclusion  of  this  section,  it  is  clear  that  modelling  usage  by  a  real-valued 
stochastic  process,  instead  of  a  Poisson  process,  may  make  matters  more  realistic 
but  does  increase  the  complexity  of  the  resulting  density.  The  situation  examined 
in  this  section  is  a  good  example  of  what  can  happen;  the  compound  Poisson 
process  is  one  of  the  more  analytically  tractable  real-valued  stochastic  processes 
vet  it  is  still  difficult  to  define  the  associated  joint  density  of  failure. 


Figure  4  Simulated  joint  density  for  (4.7) 


5  The  usage  process 

As  an  aside  to  the  probability  model  proposed  in  this  thesis,  an  idealised 
stochastic  process  for  modelling  usage  of  an  automobile,  computer  or  any  device 
where  usage  is  a  real-valued  quantity  is  introduced.  This  process  will  be  called 
the  usage  process.  The  usage  process  is  a  function  of  three  sets  of  non-negative 

valued  random  variables  : 

i)  Xv  X2,  X3,  ...  define  the  length  of  consecutive  periods  of  rest  and 

ii)  Y1,  Y2,  Yg,  .-  define  the  length  of  consecutive  periods  of  use  and 

iii)  Gv  G2,  G3,  define  the  rate  of  use  during  each  period  of  use  Y2  etc. 

The  process  starts  with  rest  period  Xr  At  the  end  of  this  period  the  first  usage 
period  Y-.  starts,  with  usage  accumulating  at  the  random  rate  G^  When  Y-^ 
finishes,  X2  begins  and  so  on,  with  each  rest  period  Xj  followed  by  a  use  penod 
Y-  During  the  ith  usage  period,  amount  of  use  increases  at  a  rate  Gj.  During  rest 
periods  the  amount  of  usage  stays  constant.  The  amount  of  use  in  the  i-th  use 
period  will  be  Y^.  Figure  5  shows  a  typical  realisation  of  the  usage  process. 

Let  M(t)  denote  the  amount  of  use  at  time  t.  First  of  all,  if  we  suppose  that  at 
time  t  we  are  in  the  (n+l)th  rest  period,  then  usage  to  time  t  is  given  by 


(5.1) 


M(t)  =  £  YjGj 
i=l 


Now  suppose  at  time  t  we  are  in  the  (n+l)th  use  period,  which  started  at  time 

T  ,  Then  the  usage  to  time  t  is  given  by  the  following  equation  : 
n+1' 


M(t)  =  YiGi  +  ^  '  Tn+l)Gn+l 
i=l 

where  T^  =  X^ 

Ti  =  Ti-i  +  yh  +  xi 


(5.2) 


i— 2,  3,  ... 
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We  can  combine  (5.1)  and  (5.2)  into  one  expression  by 


(5.3) 


M(t)  =  52  ^i’  ^(t'Ti)  x(t‘Ti)  } 

i=l 

where  are  defined  as  in  (5.2) 


and  h(x)  = 


{ 


1 ,  if  x  >  0 
0,  otherwise 


Figure  5  Typical  realisation  of  the  usage  process  M(t). 

In  practice,  one  could  allow  one  or  more  of  the  three  sets  of  random  variables  to 
be  deterministic  if  the  situation  required  it;  i.e  if  M(t)  measured  the  length  of 
time  a  light  bulb  is  switched  on  up  to  time  t  then  G-l  Vi  would  be  sensible.  If 
the  X-  are  independent  and  identically  distributed,  and  the  Yj  are  too,  then  the 
sequence  (Xp  Yp  X2,  Y2,  ...)  is  an  alternating  renewal  process. 
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6  Conclusion 


This  thesis  is  an  initial  attempt  at  addressing  some  of  the  requirements  for 
models  indexed  by  time  and  usage.  The  approach  has  been  to  make  assumptions 
about  the  relationship  between  the  two  scales  of  measurement  and  to  develop  a 
probability  model  from  those  assumptions.  In  section  2,  as  well  as  obtaining  a 
class  of  joint  densities,  a  class  of  marginal  densities  for  time  to  failure  was 
obtained  that  may  be  considered  a  viable  univariate  model  for  failure  in  its  own 
right.  By  allowing  specification  of  a  rate  of  use  function  A(t)  and  a  hazard 
function  r(t)  to  describe  environmental  effects,  the  model  offers  great  flexibility. 
The  discrete  nature  of  the  usage  variable  U  is  not,  I  believe,  a  cause  for  concern 
since  one  is  at  liberty  to  make  the  unit  of  use  relatively  small  (so  for  example  in 
the  case  of  automobile  mileage,  setting  U  to  be  the  number  of  miles  driven). 

Section  4  considered  the  replacement  of  the  Poisson  process  to  model  usage  by 
some  other  stochastic  process,  and  gave  an  example  which  showed  that  one  can 
expect  greater  problems  in  deriving  the  model.  Generally,  my  attempts  to  make 
usage  a  continuous  quantity  have  met  with  only  partial  success  so  fax. 

Section  5  considered  a  more  realistic  stochastic  process  to  model  usage  than  is 
used  in  the  rest  of  the  thesis,  named  the  usage  process.  Ideally  one  would  like  to 
obtain  a  model  of  the  form  proposed  here  that  used  this  process  for  M(t),  but 
there  remains  a  great  deal  of  investigation  to  be  done  before  that  can  be 
achieved.  It  is  quite  possible  that  analytic  results  will  be  unobtainable,  although 
one  should  be  able  to  simulate  the  process  if  this  is  the  case.  It  is  my  intention  to 
continue  research  on  this  particular  topic. 

In  summary,  then,  this  is  a  new  area  of  research  and  there  are  many  avenues  of 
thought  still  to  be  explored.  There  are  entirely  different  methods  of  creating 
models  indexed  by  two  variables  that  need  to  be  investigated,  and  the  method 
used  here  needs  to  be  extended  to  other  usage  processes,  perhaps  to  other  forms 
for  the  conditional  failure  rate  and  even  to  looking  at  models  indexed  by  three  or 
more  dependent  variables.  In  the  immediate  future,  my  research  in  this  area  will 
concentrate  on  the  usage  process  and  on  trying  to  develop  a  model  around  it. 
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Appendix 

Deriving  a  failure  model  indexed  by  time  and  usage 

This  appendix  shows  how  the  failure  model  indexed  by  time  and  usage  may  be 
derived.  It  will  consider  the  specific  cases  that  were  introduced  in  sections  2  and 

4. 

Recall  from  equation  4.3  that  we  have  a  general  expression  for  the  density  of 
failure  at  time  t  and  usage  u,  when  the  conditional  failure  rate  is  r(t)+/?M(t)  and 
M(t)  is  a  stochastic  process  describing  usage. 

t 

(Al)  f(t,  u)  =  f(M(t)=u)  x  (r(t)+/?u)  e'R^  exp  (-/?  j  M(s)  ds  |  M(t)=u  )  . 

0 

In  the  following  two  sections  we  will  show  the  derivation  of  (Al)  when  M(t)  is  a 
Poisson  process  (as  in  section  2)  and  a  compound  homogeneous  Poisson  process 

(as  in  section  4). 

The  Poisson  process 

Let  M(t)  be  a  Poisson  process  with  intensity  function  A(t).  Its  density  is  Poisson 
with  mean  A(t).  What  about  the  integral  that  equation  (Al)  says  we  must 
specify  to  find  the  density?  Let  Sj  denote  the  time  of  the  i-th  event  in  the 
process.  If  M(t)=u  then  the  area  under  the  process  from  0  to  t  is  : 

t  u 

(  f M(s)  ds  |  M(t)=u  )  =  tu  -  J2Si 

0  1=1 

Conditional  on  M(t)=u  then,  as  an  unordered  set,  standard  Poisson  process 
theory  tells  us  that  the  Sj’s  are  distributed  as  a  random  variable  S  which  has 
density  A(s)/A(t)  on  [0,  t),  and  is  0  otherwise.  So 


=  E  (  exp  (-/?( tu-Su))  |  M(t)— uj 
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u 


=  |  E  (  exp(-/?(t-S))  I  M(t)=u  )  ]•  u 

|  exp(-/?(t-s))  .  A(s)/A(t)  ds  | 

0 

Substitute  this  expression  into  (Al)  to  obtain  equation  (2.4)  as  required. 


Compound  Poisson  Process 

In  section  4  we  attempted  to  generalise  the  model  of  section  2.  To  simplify 
matters  a  little,  we  assumed  that  r(t),  the  non-use  failure  rate,  was  a  constant  a. 

Let  M(t)  be  a  homogeneous  compound  Poisson  process  with  intensity  A,  and  let 
each  jump  in  the  process  be  an  independent  and  identically  distributed 
exponential  random  variable  with  mean  1.  Let  N(t)  be  the  number  of  events  to 
time  t  (so  that  N(t)  is  Poisson  with  mean  At).  Then  M(t)  can  be  defined  as 

0,  if  N(t)=0 

N(t) 

E  Uj,ifN(t)>l 
i=l 

where  the  lb  are  iid  exponential  random  variables  with  mean  1. 

The  density  for  M(t)  can  be  calculated  by  conditioning  on  N(t);  if  N(t)=n  then 
M(t)  is  the  sum  of  n  iid  exponentials  with  mean  1,  so  has  a  gamma  distribution 
with  shape  n  and  scale  1.  M(t)  can  take  any  value  in  [0,  oo).  First  note  that 

P(M(t)=0)  =  P(N(t)=0)  =  e'At 

so  that  there  is  mass  at  M(t)=0.  For  M(t)=u>0  there  is  a  continuous  density 
function 

(A3)  f(M(t)=u)  =  £  f(M(t)=u  |  N(t)=n)  P(N(t)=n) 

n=l 
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oo  un-le-u  (At)a  At 

(n-1)!  '  n! 
n=l  v 


=  e-u'At 


OO 

E 

n=l 


nn-l(Xt)n 
(n-1)!  n! 


So  for  this  paxticulax  compound  Poisson  process  we  have  a  density  function. 


As  well  as  the  density  of  M(t),  we  must  specify  an  exponent  of  the  integral  of 

M(t).  This  was  easy  when  M(t)  was  a  Poisson  process,  but  now  this  is  more 

difficult.  Let  Si  be  the  epoch  of  the  i-th  event,  and  let  lb  be  the  amount  of  usage 

at  the  i-th  event.  By  observation  of  figure  3,  one  can  see  that  the  area  under  the 

process,  conditional  on  N(t)=n,  is 

t  n 

(A4)  (  M(s)  ds  |  M(t)=u,  N(t)=n)  =  £  (t-Sj)  Uj 

L  i=l 


In  the  Poisson  process  case  M(t)=N(t)  but  here  that  is  not  so.  If  M(t)-0  then 
obviously  JM(s)  ds=0  so  that  f(t|  M(t)=0)  becomes  ae'at.  For  M(t)=u>0 
t  ,  N(t) 

(A5)  exp  {-/?{M(S)  ds  |  M(t)=u}  =  E  (exp  {  -/?  .£  (t  -  Sj)  U-  }  |  M(t)=u). 

1-1 

Condition  on  N(t)  and  average  with  respect  to  the  distribution  of  N(t)  given 

M(t) 


=  £  E  (  exp{  •/»£(*  -  Si)  Ui}  I  N(t)=n)  x  P(N(t)=n  I  M(t)=u)  . 

n=l  i=l 

Conditional  on  N(t)=n,  then  M(t)  =  E  Uj  =  u,  so  we  can  say: 

i=l 


=  £  E  (  exp{-  /?t£ui  +  |  M(t)=u,  N(t)=n)  x  P(N(t)=n  |  M(t)-u) 

n=l  i=l  i=l 

=  £  E(exp{-/?tu  +  ^£siUi}|M(t)=u,N(t)=n)  x  P(N(t)=n  |  M(t)=u) 
n=l  i=l 

=  e"^u  £  E  (  exp  {/?£siUi  }  |  M(t)=u,  N(t)=n  )  x  P(N(t)=n  |  M(t)=u) 
n=l  i=l 


Equation  (A6) 
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What  is  P(N(t)=n|  M(t)=u),  that  appears  in  (A6)?  If  u=0  then  N(t)-0  with 
probability  1,  so  P(N(t)=0  |  M(t)=0)  =  1. 


If  u>0  then  we  have  the  following  : 
(A7)  P(N(t)=n|  M(t)=u) 


P(N(t)=n)  P(M(t)=u  |  N(t)=n) 
P(M(t)=u) 


(At)D  At 
n! 


-At-u 


00 

E 

i=l 


u  e 

•  (n-1)! 

u^At)1 

(i-1)!  « 


uI1'1(At)n 

(n-1)!  n! 

^  u1'1  (At)1 


n=l,  2,  3,  .... 


It  can  be  shown,  in  a  similar  fashion  to  (A7),  that  conditional  on  M(t)-u>0  and 
N(t)=n  the  joint  density  of  (U^,  ...,  Un)  is 


(A8)  f(Ul>  nn  |  M(t)=u,  N(t)=n)  =  ^  i  f°r  >  0,  =  u 

and  that  the  joint  density  of  (S^,  ...,  Sn)  is 

(A9)  f(sp  sn  |  M(t)=u,  N(t)=n)  =  5  for  0  <  s1  <  s2  <  ...  <  sn  <  t 


and  that  the  two  sets  of  random  variables  (Ulv..,Un)  and  (Sj,...,Sn)  are 
independent.  Use  (A7),  (A8)  and  (A9)  in  equation  (A6)  to  say 


exp  |-/?|m(s)  ds  |  M(t)=u| 

0 

unl(At)n 

-Btuf*  (n-1)!  n!  r  (n-1)!  n! 

“  C  nt'i  oo  Ku^t11 

Si  (i-1)!  i! 


exp  {/?£  snj  d(s,  u) 
i=l 
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(A10) 

e'^tu 

CO 

*  £ 

oo  u1'1^)1 

ii'l  (i-D!  i’- 

n=l 

where  f2n  = 

{(sp  •••,  sn>  u]p 

...,  un)  G  R' 

Use  (A10)  and  (A3)  in  (Al)  to  give  the  joint  density  expression  of  (4.4)  that  we 
require. 


Simulation 

Numerical  integration  of  (A10)  is  a  forbidding  task.  Simulation  provides  a 
quicker  method  of  approximating  (A10),  via  equations  (A5),  (A7),  (A8)  and 
(A9).  A  simulation  strategy  to  calculate  the  density  at  one  point  (t,  u)  is 
outlined  below.  Assume  a,  0  and  A  known. 

Before  the  strategy  is  given,  there  is  one  observation  to  make.  The  distribution  of 
(U1,...,Un)  given  M(t)  and  N(t),  (A7),  is  that  of  the  order  statistics  of  (n-1) 
independent  uniform  random  variables  on  [0,u].  It  is  easy  to  show  that  if  we  take 
the  order  statistics  of  (n-1)  samples  from  a  uniform  distribution  on  [0,u],  then  the 
distance  between  each  consecutive  order  statistic  gives  a  sample  u1,...,un. 

Select  (t,  u),  where  u>0  (we  know  the  joint  distribution  for  u=0).  The 
expectation  (A5),  which  forms  part  of  the  joint  density,  is  estimated  by  the 
mean  of  a  sample  of  size  K,  where  K  is  some  suitably  large  number. 


Repeat  :  k=l 

-  using  (A7),  sample  from  the  distribution  of  N(t)  |  M(t)=u  to  find  a  value  n  . 


-  using  (A8),  sample  from  (Up  —jU^*)  I  M(t)— u, 
samples  from  a  uniform  on  [0,  u],  say  (vp...,vn*  ^ 


v(l)<v(2)<-<V-l)' 

required  sample. 


Now  u1=v^,  u2=v^-v^, 


N(t)=n*.  Generate  n*-l 

).  Order  them,  so  that 

...,  u  *=u— v,  *  as  the 
n  (n  -1) 
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-  using  (A9),  sample  from  (S1,...,Sq*)  |  M(t)=u,  N(t)=n*.  Generate  n*  uniform 
samples  on  [0,t].  The  ordered  sample  is  s^  s^*. 

* 

-ek=  expj^EsjUj}. 
i=l 

-  k=k+l.  Until  k=K 


Using  (A3)  in  (Al)  we  can  write  the  joint  density  of  T  and  U  as 

-1  t 

f(t,  u)  =  e'u'At  £  (^;)(i1U1Vx(q^u)  e'atE  (exP  {  ’  M(S)  ds)  1  M(t)=U  ) 


n=l 


0 


=  (a+/?u)  e 


-u-(a+A)t 


oo 

E 

n=l 


(At)p  u11'1 
n!  (n-1)! 


N(t) 

x  E  ^exp  {  -0tU+0 £  SjU;}  |  M(t)=u). 


Now  use  the  mean  of  elv..,ejj  to  estimate  the  expectation  that  occurs  in  the 
above  expression: 


(All)  f(t,  u)  a  (oWu)  e-(«+J)‘  £ 

n=l 


x  e 


,-/?tu  J_  ^ 


K 


1=1 


All  the  three  sets  of  samples  to  be  made  -  n*,  (ulv..,un*)  and  (s^...^*)  -  are 
easy;  the  first  is  a  discrete  distribution  and  the  latter  two  only  require  sampling 
from  a  uniform  distribution  followed  by  some  sorting.  Figure  5  gave  a  plot  of  a 
simulated  density  when  «=/?=A=l  using  the  above  procedure.  In  this  case 
P(U=0)=0.5,  so  the  volume  under  the  jointly  continuous  part  of  the  density  is 
only  0.5.  For  this  simulation  K=10000  and  the  infinite  sum  present  in  (All)  was 
truncated  at  n=100  (the  sum  converges  very  rapidly).  Calculations  on  a  PC  took 
around  4  seconds  per  density  value. 


26 


References 


Cox.  D.R.  (1972)  Regression  models  and  life  tables  (with  discussion).  J.  R. 
Stat.  Soc.  B,  v.  39,  pp.  86-94. 

Farewell,  V.T.  and  Cox,  D.R.  (1979)  A  Note  on  Multiple  Time  Scales  in 
Survival  Analysis.  Appl.  Statist.,  v.  28,  no.  1,  pp.  73-75. 

Hougaard,  P.  (1987)  Modelling  Multivariate  Survival.  Scand.  J.  Statist.,  v.  14, 
pp.  291-304. 

Lemoine,  A.J.  and  Wenocur,  M.L.  (1986)  A  Note  on  Shot-Noise  and  Reliability 
Modeling.  Operations  Research,  v.  34,  no.  2,  pp.  320-323. 

Oakes,  D.  (1988)  An  Equivalency  Model  for  Multiple  Time  Scales  in  Survival 
Analysis.  Technical  Report  88/04,  Dept,  of  Statistics  and  Division  of 
Biostatistics,  University  of  Rochester,  Rochester,  NY. 

Singpurwalla,  N.D.  and  Youngren,  M.A.  (1990)  Failure  Models  under  Dynamic 
Environments.  Technical  Report  GWU/IRRA/Serial  TR-90/2,  The  Institute  for 
Reliability  and  Risk  Analysis,  The  George  Washington  University. 


27 


